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Introduction 


In the design of a structure it is important to know how the 
structure will respond to forces that are expected to be present in its 
operating environment. Much time is spent developing analytical models 
that predict responses, and when the structure is finally built it needs to 
be physically tested before being placed into operation, This modal 
testing has become very important in the last few decades with the struc- 
tures becoming more complex and being placed in more extreme operating 
environments. 

In Chapter 1 of this thesis, a review of several modal testing 
techniques will be made, along with brief discussions of their advantages 
and limitations. In Chapter 2, a new technique will be presented which 
overcomes many of the previous limitations. Several simulated experiments 
are Included in Chapter 3 to verify the validity and accuracy of the new 
method. Finally, in Chapter 4 conclusions are drawn from the simulation 
studies and recommendations for further work are presented. 

The complete computer code configured for the simulation study is 
presented in an Appendix. 
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Chapter 1 


1.1 Overview 

The need for modal testing (sometimes referred to as modal analy- 
sis) of structures is well known and documented. The problem is that there 
exists no fool-proof and totally efficient method to extract modal coeffi- 
cients from structures. Over the past 40 years several different types of 
analysis have been developed and practiced with varying degrees of success. 
The analyses range from time-domain to frequency-domain, and from single- 
excitation to multiple-excitation. To employ each of these methods proper- 
ly, it is necessary to examine them to find their respective advantages and 
disadvantages. 

1.2 Multi -Shaker Sine Dwell Testing 

The oldest and most widely used method is the Multi-Shaker Sine 
Dwell (MSSD) test, first described by Lewis and Wrisley (1). This method 
uses shakers attached at several points on the structure. The shakers are 
tuned in amplitude and phase simultaneously at a given frequency such that 
the acceleration outputs are 90 degrees out of phase with the corresponding 
force inputs (this is indicated by a circle on a Lissajous plot). The 
amplitudes and phases are then recorded giving the mode shape of the 
structure at that given frequency. The damping can also be obtained by 
simultaneously shutting off all the applied forces and analyzing the 
resultant free-decay rate. 
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The advantages of this method are that a very accurate mode shape 
can be obtained due to the high signal-to-noise ratio, and Input levels can 
be varied to permit Identification of pseudo-modal parameters for non- 
linear structures. The disadvantages, though, usually outweigh these 
benefits. To measure the mode shapes accurately, accessibility to several 
shakers Is required, and In some Instances the attachment of a shaker at a 
given location is physically impossible, A second drawback is the time 
involved In stepping through the frequency range one mode at a time and 
tuning and re-tuning each shaker to obtain a 90 degree phase lag (or lead), 

A revision of the MSSD method by Anderson (?.) speeds up the tuning 
procedure. In this method, interest is placed on maximizing the quadrature 
part of the output to input relationship. Several closely-spaced modes are 
tuned in conjunction with one another with the number of shakers used being 
dependent upon the number of mode? being tuned at a tv-*,. Though the time 
required by this method is less, it is still an iterative procedure, which, 
in itself, makes this method time consuming. 

1.3 Time-Domain Analysis 

As far as time-domain analyses go, there are basically two types, 
the Least-Squares Curve-Fitting Time-Domain (LSTD) Method (3) and the 
eigenvalue Ibrahim Time Domain (ITD) Method (4). The LSTD method is an 
older and more widely known method than ITD, though ITD is more sophis- 
ticated and has fewer limitations. 

The LSTD method uses the time history of a free-decay response and 
upon an initial approximation of the modal parameters, it performs an 
iterative solution routine, based on the Newton-Raphson procedure, to 
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produce accurate estimates. The problems with this method are the computer 
time Involved for the Iterations {up to 50), the need for user-supplied 
Initial approximations, and the reported high Inaccuracies after the first 
two or three modes (3). 

A revision of the L$TD method has recently been developed by Smith 
(3), In this revision a new type of Newton-Raphson iterative procedure is 
used. This method Is also semi-automatic, with the computer analyzing an 
FFT of the free-decay response to gather the initial natural frequencies 
(the damping estimates are still user-supplied). This revision, though, 
also has difficulty in separating closely-coupled modes. If two peaks 
cannot be resolved in the initial FFT, then one of the two modes will be 
omitted from the analysis. Also, to properly Identify two closely-spaced 
modes, the time history must Include one or two oeat frequency cycles, thus 

i 

causing a significant increase in the length of the required data record. 

The ITD method uses time histories of free-decay responses and 
time-delayed free-decay responses of several degrees-of-freedom (DOF) to 
set up an eigenvalue problem that solves for the modal coefficients. To 
solve for more modes than there are measured DOF, ITD creates pseudo-DOF by 
using time-delayed responses based on the original responses. This also 
provides the means for determining the validity of each acquired mode by a 
comparison of morie shape coefficients. 

The ITD method has proven to be an effective tool in modal analy- 
sis. There is a problem, though, in the use of this method: it is fairly 
complicated. For each test there are seven analysis parameters to choose. 
Also, when there is a low signal -to-noise ratio, ITD produces some system- 
atic errors that must be compensated for by several additional analyses 
(5). 
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1.4 Frequency-Domain Analysis 

Frequency-domain analysis encompasses several different styles of 

curve-fitting or eigenvalue analysis. The common features among all the 

methods are their similarity of excitation and their use of curve-fitting 

« 

of frequency response functions (FRFs). The structure Is commonly excited 
by one shaker with an input signal in the form of a swept sine wave or some 
type of random noise. Sometimes, though, the structure may be excited by a 
single impulse, usually administered by a hammer with a force cell at- 
tached. 

FRFs are usually acquired by a single-input, single-output proce- 
dure. The system is either excited by a hammer or a single point random 
(SPR) shaker set-up (6). By dividing the cross-spectrum of the input and 
output by the auto-spectrum of the input an accurate estimate of the FRF 
can be obtained. 

Another procedure for gathering transfer functions has been de- 
scribed by Allemang (7). In this method the structure is excited simulta- 
neously by several random-noise-driven shakers. The excitations and 
responses are measured and then manipulated. 

{x} = [HI {F} 

[GXF] = [H] l GFF] 

[H] = [GXF] [GFF]’ 1 


where, 

[GXF] = m x n frequency dependent matrix 
[GXF] = Input-output cross-spectrum matrix 
[GFF] = n x n frequency dependent matrix 
[GFF] = Input cross-spectrum matrix 

GFF“( i » J ) = Cross-spectrum of input i with respect to input j 


C 
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This method has the advantage over SPR In that several columns of 
the transfer matrix can be computed at one time. This method has the 
disadvantage, though, of a high computational ’cost,' since the [GFF] 
matrix must be Inverted at each frequency. 

The simplest curve-fitting technique is the 'peak-picking* method, 

A plot of Imaginary response vs, frequency Is used to find the points of 
maximum magnitude in the Imaginary component of the response. If the 
structure can be represented as a single DOF (SDOF) system In the vicinity 
of the natural frequency, the point of maximum magnitude corresponds to the 
natural frequency, and its value is the mode shape component. A half-power 
scheme can be employed to determine the damping value. 

If another mode is close to the chosen mode, errors in the modal 
parameters will become significant. One way to compensate for this is the 
SDOF circle-fit technique. With this method a plot of imaginary vs, real 
values of the FRF is analyzed. When a natural frequency occurs, a circle 
will tw formed in the Argand plane. If two or more modes interact, the 
circles will not be complete; one circle will be only partially formed when 
the next circle starts (8). An estimate for each complete circle can be 
made, and from this the amplitude of the mode shape (diameter of the 
circle), natural frequency (mid-point on the circle), and the damping (the 
circle is bisected and the corresponding magnitudes and frequencies are 
used) can be obtained. This procedure was originally suggested by Kennedy 
and Pancu (9), 

This SDOF circle-fitting method will also break down if the modes 
become too closely spaced such that the individual circles interfere 
significantly with each other. 


ft 


n 

One method that can help produce better estimates is a routine 
described by Richardson and Kniskern (10), With this method FRFs are taken 
from the structure using various reference exciter locations, so that 
components from more than one column of the transfer matrix, H, are in- 
cluded, After initial analysis, the natural frequencies and residues are 
input to an algorithm that weights and averages the residues to produce 
modal vectors. 

The problem with this method is that it is not self-sufficient nor 
analytically rigorous. Since a preliminary analysis must be performed 
beforehand, the time factor Increases significantly. The added accuracy is 
not always sufficient to warrant the use of this algorithm either. A 
numerical example cited in the paper by Richardson and Kniskern showed only 
a slight Improvement as a result of using this method, 

There are, though, several curve-fitting routines that can properly 
handle multiple DOF (MDOF) situations. Among them is one called the 
Least-Squares Complex Exponential (LSCE) technique (11). In this method 
the FRFs are obtained and then transformed back into the time domain, 
forming the unit Impulse response functions. The analytical form of the 
response function can be written 


X(t k ) 


2N 

£ 

r=l 



where, 


A r - residue of mode r 

A r+ i * complex conjugate of A r (r=l ,3,5,. , . ) 



(l-i) 


ORIGINAL PAGE ?3 
OF POOR QUALITY 

12 

s r * -VV * 

s r+1 * complex conjugate of s r (r*l,3,5,..,) 

T * value of time sub Interval 

k * time point index (k-0,1,2,. * „ ,2n) 

2n * number of time steps 

N * number of modes 

If N ■ n, then there is a polynomial of order 2N that must 
satisfy the equation 

U 2N + a^" 1 + a^” 1 + ... + a^Uj + a 2N * ° (1-2) 


whose roots are Up U 2 ,.,., , 

Pre-multiplying Eq, (1) by §n - summ * n 9 ytelds 


2N 


2N 


2N 


? n a 2N-k X(t k> = ^ A r { ,L a 2N-k U r > 


k~0 

The following equality holds 

2N k 

k ? 0 a ZN-k U r = 0 


k®o 


since U r is a root of Eq. (1-2). 
Thus 


2N 

k s0 a 2N- k “ ° * a 0 s 1 

Equation (1-3) can be expanded into a series of equations by beginning each 
successive equation with a different response time, i.e. 
m 2N 

1=0 ^k=o a 2N-k x ^k+l' = °) * m >2N 0 - 4 ) 

Since a Q is known, the other a’s may be obtained from the equation 




* 


•«w> 
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(1-5) 



2N-1 

^k*0 a 2N-k X ^t+1 


Once the a's are known they can be substituted back Into Eq. (1-2) and a 
polynomial root solver can be used to find the roots, U M , The values for 
w r and can be found directly from U r , The residues, A r , for 
each Up can be found by creating 2N equations In the form of Eq. (1-1), 
each with a different Initial time* 

Though this method was just described for one unit impulse response 
record, It can be easily extended to several records. Since the a's are 
system parameters, they should remain constant for each response record. 
Therefore, the summation on i can be extended to Include any number of 
response records. This also means that the excitation location can be 
changed so that more than one column of the H matrix can be used simulta- 
neously In t!)f TOlysis* 

To a^Labllsh whether the estimate for N , the number of modes, is 
sufficiently large, a simple procedure is used Involving a series of 
progressively larger estimates for the number of modes. Equation (1-5) is 
first rewritten in matrix form 


Aa = B { 1-6) 

and the least-squares error Is found by computing 

R B Aa - B (1-7) 

and the sum square error E Is 
= r t r 


E 


(1-8) 


A plot of the error vs. the number of modes assumed will show a significant 
drop-off when N Is sufficiently large. 

The LSCE method has been incorporated Into commercial software 
packages and has demonstrated great success In extracting modal parameters. 
LSCE does have a problem, though. It must transform frequency-domain data 
back to the time domain. When this additional transform is performed, 
Information can be lost or contaminated. 

A second method which uses FRFs from SPR testing and which can use 
multiple columns of the transfer matrix, H , is an eigenvalue solver 
developed by Coppolino, called the Simultaneous Frequency Domain (SFD) 
Method (12). The idea of 'simultaneous' is that single-input FRFs from several 
response locations and reference locations can be used simultaneously in 
the analysis. In this method the matrix equation of motion for a structure 
may be transformed to the frequency domain and written as 

CM] { X ( w ) } + [B] {X(w)} + [K] (X(w)} = {D} f(w) (1-9) 


where 

[M] = mass matrix 

[B] s damping matrix 

[K] - stiffness matrix 

(X(w)} = array of displacements of points on the structure 
{D} - external force distribution array 

f(w) = external force amplitude 

f(w) - 1, since the transfer function for a unit impulse is 1. 
Pre-multlplying Eq. (1-9) by [M] and solving for (X(w)} gives 
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% 


{X(w)> 


-[CM ]' 1 [B] | [M ]' 1 [K] | CM ]- 1 CD]] 



( 1 - 10 ) 


The response arrays can be expanded Into matrices such that each 
column corresponds to a frequency point. A least-squares solution can then 
be obtained for the quantities [M]" 1 [B] and [M]" 1 [K]. Once these two 
dynamic matrices are known, it Is a simple procedure to use them in an 
eigenvalue problem 


fx' 

' + 

" Cm] -1 [b] 

[M]' 1 [K] ' 


r • 

X 

► =5 . 

' o' 

l*J 


. -D3 

CO] 


X 

* * 


0 


to solve for the frequencies and modes. 

The only disadvantage of SFD compared to the LSCE single-input MDOF 
parameter solver Is that the order of the matrices, thus the number of 
modes obtained at a time, is directly dependent upon how many response 
locations are available for the analysis. 

A third method that utilizes multiple rows and columns is the force 
appropriation, or Asher, method (13-15). The Asher method is based on the 
same idea as the MSSD test; that at a natural frequency an appropriate 
force vector will produce responses 90 degrees out of phase with the 
excitation. Assuming that the applied forces are real, this means 


{X} - [H] {F> = ([H r ] + 1[H,]){F> (1-12) 

For the response to be in quadrature with the excitation 

[H r ] {F} = 0 (1-13) 


Thus, the condition for natural frequencies is 
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Det [H r (w r )] = 0 r « 1,2 n (1-14) 

Once the frequencies have been found from Eq. (1-14), a solution 
for the required force vector, J*“ r , can be obtained by using Eq. (1-13). 
Finally, the mode shape can be obtained from 

{X(w r )}= i[Hj (vv r )]{F> (1-15) 

The Asher method is usually accurate, but the condition of Eq. 
(1-13) is only a necessary, but not sufficient, requirement for a natural 
mode. If the order of the determinant is much less than the order of the 
system, then artificial or 'spurious* modes can appear. 

The Asher method requires that the number of excitation and re- 
sponse points be equal, whereas the number of response points usually far 
exceeds the number of excitation points used in a test. To avoid this 
Asher method limitation, another method has been developed which usee 
non-square matrices, and thus permits more response points. This method, 
termed the minimum coincident response (MCR) method (15-16), starts with 
basically the same equations 

{X} ~ ([H R ] + 1[Hj]) {F} p«n m<n (1-16) 

m x 1 m x p p x 1 

Or} » [Hr] {F} (1-17) 

Instead of setting Eq. (1-17) to zero, a coincident response error norm 
given by 

e = M T {M 

is minimized. The force vector F is appropriated such that e is 
minimized subject to an amplitude constraint on j , ' ** quadrature 
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response at the ith response point. This leads to a force vector 


« - 


(M Mr? CV)' 1 ( h ii) 


([H r 3 T [H,])- 1 {H n } T (1-18) 


where Hjt is the ith row of the imaginary part of the transfer matrix, 
The mode shape can then be obtained from 


{X,} = i[H i (w r )] {f} 


(1-19) 


The MCR method has been shown to give favorable results (15), 
although it can still produce 'spurious' modes. This problem can be 
diagnosed by various methods. One is to check the slope of the e vs, 
frequency curve. Shallow slopes seem to indicate 'spurious' modes, while 
steep slopes correlate with actual modes (15). Another difficulty with the 
MCR method is that, while the Asher method interpolates to find zero 
crossings, MCR searches for minima. 



Chapter 2 


2,1 Development of a New Multi-Shaker Algorithm 

For several years, since the advent of digital computers and the 
development of the Fast Fourier Transform (FFT), considerable time has gone 
into performing and perfecting single point random (SPR) excitation struc- 
tural dynamics measurement techniques. The SPR method has proven to be 
reasonably fast and simple to perform. There are, though, some inherent 
problems with an SPR test. Two of the chief problems are energy dissipa- 
tion within the structure and the possibility of locating a shaker near a 
node line. 

A technique to avoid these problems is to apply forces at multiple 
locations. The method to be developed in this chapter is an expansion of 
the eigenvalue SFD method described by Coppolino (12). Since the standard 
modal analysis methods are applicable to single-shaker testing only, this 
alternative approach involves a reformulation of the equations of motion to 
accommodate the use of multiple forces. 

The Multi-Shaker Modal Analysis (MSMA) method presented here 
assumes that a linear structure can be accurately described over a limited 
frequency range by a finite number of degrees of freedom (DOF). This 
assumption yields the standard frequency-domain matrix equation relating 
the output acceleraticns at given points on the structure to the input 
forces 

{X(w)> = [H(w)] {F (w )} (2-1) 
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where, 

f «• 

\X(w)j = output accelerations , (n x 1) 

(f(w)} *» input forces , (n x 1) 

CH(w)] = complex frequency response function , (n x n) 

[H (w ) ] - [0] [/A/] [0] T 

[*0] a modal matrix , (n x n) 

[/A/] = diagonal modal admittance matrix , (n x n) 

C/A/] = (Iw/-- + f(2c r «w r )) r=i ’ n 

W r = natural circular frequency of mode r 
£ r = critical damping ratio of mode r 
w = circular frequency 


The force vector is rewritten in a form that includes only the DOF 
where forces are applied. A distribution matrix relates this new force 
vector, (f{w)} , to the original force vector, {F(w)} . 

{f(w)} = [D] (f(w)} (2-2) 


For example, if a 4D0F system has exciters at points 1 and 3, the equation 
would look like 



f V 

o 1 


1 

0 


f 

= 

0 

0 


k 

[ f 3j 

fo 


0 

1 


0 3 J 


0 

0 



So, Eq. (?-l) is now 

{X(w)} = [0] [/A/] [0] T [D] {f (w)} 


(2-3) 
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It is assumed that there are n points where the responses are measured, 
and if it is assumed that the structure can be accurately described over 
the frequency range of interest by m DOF (m < n), then Eq. (2-3) can be 
partitioned into sets of 'independent' DOF and 'dependent' DOF equations. 
Thus, 


5c,(w) ’ 


C0,]' 

r t t1 

[o,] 

= 

C/A(w)/3 {[0 f ] 7 [0 d ] 7 J 

[k d (w) j 


_[0 d ] _ 

. fD d ] . 


{f(w)J (2-4) 


or 


'*>) ' 


‘ m 

V*> 

J 


CP d ] [ 0,]* 1 


CP,] [MM/] [0/ [D 1d ] {f(w)} (2-5) 


where 


[D 1d ] - [D-] + W,]' 7 C0 d ] T [D d ] 

From Eq. (2-5) a relationship between the independent and dependent DOF 
accelerations is 

{x d M} = [p d ] [0 ,]" 1 {)(,(«)} (2-6* ) 

This equation will be employed later to determine the modal components of 
the dependent DOF for the determined natural frequencies. The equation 
relating the independent accelerations to the forces is 

{if,.} ■ [0,] C/A(w)/] [0,.] T [D 1d ] {f(w)j (2— 6b ) 

Pre-multlplying both sides of Eq. (2-6b) by ([0^.1 [/A(w)/1 [0^] 7 )’ 7 and 
expanding, yields 
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([^3 c/a(w)/ 3 c^rr (xj ■ [D 1d ] (f(w)} 

l^T r [/( 1 /A(w))/J tt ,]" 1 {x'J - [D 1d ] {f(w)} 

»,r T c/d - i 2 ? r (w r /w) - (w r /w) 2 )/] cp 1 r 1 (xj - [o 1d ] (f(w)} 

[|3 1 ]‘ T [U i r 1 {5i 1 }+ [0 1 ]" T C/25 r W r /] C0,] -1 M/w) {x,.} 

+ [0 i ]‘ T [/w r 2 /] W,]- 1 (-1/W 2 ) {xj - [D 1d ] (f(w)} 


or, in summary, 


[M^ {*,} + [C,] {xj + r.K.] {xj = [D 1d ] {f (w) } (2-7) 

where 

[M-j 3 » w ,]- 7 ^r 1 

[c,] = [0 1 ]" T c/2? r w r /] 

[K t ] = [0 1 ]' T C/w r 2 /] [0,]" 1 

{*, } = (-1/W) {*,} 

{x, } - (-1/W 2 ) {X.} 

Least-squares averaging may be employed to estimate the matrices 
[M.j] , [C.] , and [K^] . Therefore, the accelerations, velocities, 

displacements and forces in Eq. (2*?) are expanded into matrices, where 
each column corresponds to a discrete frequency point. That is, 


[X,] = [X i (Wj) I X^Wj) I...I x 1 (w k )3 


(2-8) 


-1 


where k > m . Multiplication of Eq. (2-7) by [M^]” yields 


- [[?,] [K k ] [B 1d ]] 


[x,] 
Cx i .l 
L Cf(w)] 


(2-9) 


where, 


[c,] 

[K,] 

[ 6 1d ] 


[M,]’ 1 [C,] 

CM ,]* 1 [K,] 

-W l C°ld 3 


ORIGINAL PAGE IS 
OF POOR QUALITY 


22 

(2-10a) 

(2-10b) 

(2-lClc) 


Since [£,] , [K,], and [8, d ] are real matrices, Eq. (2-9) can be 

rewritten for the real and Imaginary parts of the responses and Inputs. 


[ [ V 3 tflj] = Ofy f~ K 1 3 f 5 1d 3 ] 


[X, ] 
'Re 

rx, ; 

' Im 

[X, ] 
'Re 

[X, 

'Im 

Cf Re ] 

CV 


(2-11) 


Since this equation is over-determined (k > m), a least-squares solution is 
used. Elation (2-11) may be rewritten symbolically 


[Y] - [T] [Z] (2-12) 

The least-sauares solution for [T] is 

[T] = [Y] CZ] T [[Z] [Z] T ] - 1 (2-13) 

The [C^] and [K^] that were solved for in Eq, (2-13) can now be 

used to determine the eigenvalues and eigenvectors for the system in the 

following manner. From Eq. (2-9) the free vibration equation 

X, + [c,3 X, + [K,] X, = 0 (2-14) 

may be written. Since the matrices [2,] and [K,] are constant, X and 

its derivatives in Eq. (2-14) may be written in the time-domain, rather 
than in the frequency-domain as originally expressed in Eq. (2-7). Thus, 


let 
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v V M 

*1 - Xe (?-X5a) 

X, * AX 1 (2-lSb) 

^ = XX, (Z-15c) 

From Eq. (2-lSb) 

AX 1 - X. = 0 (2- 16d ) 

Thus, combining Eqs. (2-14) and (2-15) yields 


X . 

[*1 

1 

. Jf. 

’rc,) 

rK t ]' 


r • i 
x i 

C 4 

0 


w 


_c-n 

Co] _ 


, x ij 


0 

» d 


which is a standard linear eigenvalue problem. The eigenvectors will 
produce mode shape components for the ‘independent 1 DOF only. To include 
the ‘dependent 1 DOF, Eq. ( 2-6a ) is expanded to matrix form for the k 
frequency points and then manipulated, via the least-squares technique, to 
produce a relationship between the ‘dependent 1 and 1 independent* modal 
coefficients. 

[XJ - W d ] nu- 1 EX.,] , % 

T T 1 (2-18) 

[ 0 d ] - [X d ] [x,] T Ccx 1 ] Cx d ] T ] - 1 [ 0 ,] 

2.2 Forcing Functions 

As indicated by Eqs. (2-4) and (2-5), the forces may be applied at 
the dependent DOF as well as the independent DOF. This feature of the 
equations allows for forces to be applied where either accelerations could 
not be measured or where an accelerometer produced bad readings and was 
subsequently discarded. It also allows the user more freedom in 


determining which locations are to be Identified as Independent DOF. If, 
after data acquisition is completed and a general analysis done, it is 
determined that a particular DOF is not Important, its responses can be 
reduced out of further calculations, even though it was forced. 

In addition to the freedom of location of forces, there is also a 
freedom in the type of forces which can be applied, subject to one con- 
straint. From Eq. (2-13) it is seen that, if two forces are fully coherent 
the inverse of the matrix will be singular. However, the restriction is 
not too severe. If two forcing functions are Identical, only one forcing 
function needs to be included in the force array. The force distribution 
matrix will handle the duplication in the applied force spectrums. Since, 
in reality, a structure will have different input Impedances at each 
location, similar signals into the shakers will produce different forces. 
The best way to avoid any potential problems, however, is to have a low 
coherence between forces. This can be achieved by using either of two 
forcing methods, independent random signals at each shaker or simultaneous 
sine sweeps with each shaker having a different phase shift. 

Another type of excitation is a free-decay response. A minor 
reformulation of the equations is involved, with Eq. (2-11) being modified 
to the form 


DX ] EX, ]] = Cep [K-j.1] 

Re Im 


ri ] ex, ] 

’Re ’im 

EX, ] EX, ] 

’Re ’im 


(2-18) 


This formulation eliminates the need for measurement of the forces, thus 
allowing a more general type of excitation. This can also be used in 


conjunction with the Randomdec method where a free-decay response Is 
generated frorr a structure that is being randomly excited* 

2.3 Averaging 

There are two possible methods to eliminate noise contamination in 
the eigenvalue results. Both methods Involve a type of signal averaging. 
The first method Involves averaging FRFs from several runs to ff orm a single 
acceleration and force FRF spectrum at each location before insertion into 
Eq. (2-13), 

There is one serious drawback to this method. With random excita- 
tions, excessive averaging can produce an FRF spectrum with a uniform 
value, effectively giving fewer independent columns in Eq. (2-13) with 
which to do a least-squares fit. For a swept-slne test, this averaging 
technique should suffice, since the FRF will not tend to become uniform, 
especially if the signal generator has a phase modulator. 

The second method is to string the multiple runs together in the 
same manner as the individual frequencies were in Eq. (2-8), That is 

[X(w)3 - [[)(>)] | [X 2 (w)l|... | [Xp(w)J] 

This increases the number of columns, therefore the accuracy of the least- 
squares fit, but the required storage space for the matrices is also 
increased. 

2.4 Pseudo-Degrees-of-Freedom 

If the structure being tested has more modes in the frequency range 
of interest than there are measurement locations, there is a technique to 
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increase the number of DOF based on existing measurement DOF, This method 
involves sampling extra sets of data and placing each new set Into the 
equations as if they were new DOF. With this approach is an alternate 
method in which the responses of the originally measured data are 
partitioned into subsets and correspondingly labeled as different sets of 
data. In this way the number of DOF is increased without the need for more 
measurements. This, however, gives fewer columns of data for the 
least-squares process, since this procedure takes some of the columns in 
Eq. (2-8) and makes new rows out of them. 

2.5 Mode Selection 

Proper selection of the true modes of vibration from the eigen- 
solution of Eq. (2-16) is just as Important as any other part of the modal 
analysis. When the number of modes is overspecified, either through 
oversight or to compensate for noise contamination, a systematic method is 
needed to identify the spurious modes. Usually the spurious modes can be 
eliminated immediately by a quick examination of the eigenvalues. The 
damping values are often excessively high, low, or even negative. The 
natural frequencies might have values that are out of the measured frequen- 
cy range. At other times, however, the spurious mode eigenvalues are 
within the range of tolerance. Then, an inspection of the eigenvectors is 
needed. One procedure involves the use of pseudo-DOF. A comparison of the 
portion of the eigenvector representing the pseudo-DOF with the portion of 
the eigenvector which corresponds to the independent DOF from which the 
pseudo DOF were produced can readily determine the validity of the mode. A 
favorable relationship would be one where there is a high correlation of 
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eigenvector values at corresponding DOF, Another method would be to 
analyze results of analyses from several data sets; modes that repeatedly 
show up would Indicate valid modes. 

2,6 Time-Domain Analysl t 

The MSMA procedure just described may also be performed directly in 
the time domain instead of first transforming the signal Into the 
frequency-domain, although this raises some new problems, The relationship 
between acceleration, velocity, and displacement is not as simple. To 
acquire the three quantities there are three possible methods. The first 
is to measure all three quantities directly. The second is to measure the 
acceleration and then integrate to obtain the velocity and displacement. 
In this case, the initial velocity and displacement are needed. The third 
method is to measure the displacement and then differentiate the signal to 
obtain the velocity and acceleration. Since none of these procedures is 
easily accomplished, the frequency-domain approach remains the preferred 
approach. 






chapter 3 


3.1 Simulation Studies 

To verify the validity and accuracy of the MSMA technique described 
in Chanter 2, several computer simulations were developed and examined. 
The primary concerns evaluated In the tests v/ere: (1) resolution of the 
modes In high noise, (2) Identification of modes less than $ r w r apart, 
(3) accuracy of dependent DOF, and (4) application of forces at dependent 
DOF. 

A computer program was designed such that the system natural 
frequencies, damping, and mode shapes, and a nolse-to-signal ratio are 
Input. The noise level is calibrated by using the root-mean-square of the 
response signal over the entire frequency range. Also input are the 
frequency range, number of independent DOF and their corresponding lo- 
cations, number of dependent DOF and their locations, number of forced DOF 
and their locations, and the method of excitation. The program creates the 
frequency responses for the accelerations and forces and then proceeds to 
manipulate them in accordance with the equations of Chapter 2. 

The first series of tests was on a system with nine natural fre- 
quencies, seven of which were below 10 Hz. Two of the frequencies were 
selected less than £ r w r apart, so that the effects of closely-coupled 
modes could b' observed. The noise-to-signal ratio was selected at 20% in 
conjunction with an applied random force signal. Changed during the test 
series were the number of forced DOF: 1, 4, and 8 forces were used. As 
indicated in Table 1, the mode shapes had magnitudes at each station of 
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1.0, except for one DOF In each mode whose magnitude was set at 0,01 to 
simulate a point close to a node line. 

A sampling run was first made using a single shaker to obtain an 
estimate of the number of modes in the 0-10 Hz. range. A check of the FRFs 
showed that six modes were fairly apparent, so an estimate of eight modes 
(six for the observed ones and two more for any missed modes or any pseudo 
or spurious modes created by the noise) did not seem unreasonable. 

Three runs were made for each of the forced-DOF conditions plus one 
additional run utilizing the averaging technique of placing several runs 
into one large array. The averaged analysis used three individual runs 
which differed from the previous three runs. 

To check the accuracy of the results an error analysis was used 
that looked at the percentage of error in the natural frequency, damping, 
and real part of the modal vector compared to the input system characteris- 
tics. The natural frequency arid damping errors were computed directly. 
The modal vector, though, was a bit more complicated. The mode shapes all 
had magnitudes of 1 at each DOF, except for one with .01, so the estimated 
mode shapes were each scaled such that the largest component had an ampli- 
tude of 1 and a phase angle of 0 degrees. The real part of each component 
was then squared, summed (except for the one component of .01 amplitude, 
which was omitted from the error analysis), divided by 7 (the corresponding 
magnitude of the original vector); and then the square root was taken. The 
final value, which ranged from 0 to 1, was then used directly for the error 
analysis. 

To highlight the results, the errors for three modes were inves- 
tigated. The first was the 5.9 Hz. mode, which was relatively uncoupled. 


The other two modes were the coupled modes, 5,5 Hz, and 5,55 Hz. Referring 
to the diagrams in Fig. 1, It can be seen that with 20% noise the natural 
frequency estimates are highly accurate, 'aven for the closely-rrupled 
modes. Also evident in the diagrams is the overall higher accuracy 
of the parameters corresponding to the 5,9 Hz. mode compared to the 5.5 and 
5.55 Hz. modes. Although the errors for the damping and modal vectors are 
large for the 20 % noise single-shaker case, the diagrams show a steady 
trend towards lower error levels as the number of shakers is increased. In 
fact, the error levels for the closely-coupled modes are comparable to the 
error level for the 5.9 Hz. mode, when eight shakers were employed. 

Also tested with this simulation model was a set of runs in which 
the noise-to-signal ratio was changed from 0 % to 10 % and finally to 20%. 
To illustrate the effect of noise on the results and the effect of addi- 
tional shakers, the modal parameters corresponding to the 5.0, 5.5, and 
5.55 Hz. modes are tabulated (see Tables 2 through 10). An analysis of the 
tables shows, once again, a clear trend toward higher levels of accuracy at 
a given noise level for an increase in the number of shakers. A dramatic 
result from this set of runs was the absence of the 5.0 Hz mode from the 
single-shaker analysis when noise was present. As the number of shakers 
increased though, the 5.0 Hz mode appeared and then increased in accuracy. 
Overall, by a comparison of the tables, it is seen that the algorithm does 
not become unstable with increasing noise levels. On the contrary, good 
results are obtained even at relatively high noise levels when several 
shakers are used. 

The second series of tests analyzed goals (3) and (4). To do this, 
an 8 DOF system was used, with four modes in the 0-10 Hz. range. In this 
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case the noise-to-signal ratio was also 20%, but the number of forces was 
set at 1, 5, and 8 (an initial sampling run indicated three, possibly four, 
modes present in the 0-10 Hz. range, so five independent DOF were select- 
ed). 

The previous series of tests was again performed, but since the 
only substantial errors occurred in the coupled modes, their parameters 
were the only ones analyzed. In this analysis, only the normalized in- 
dependent modal vectors and a ratio of the normalized dependent to indepen- 
dent modal vectors were calculated (see Fig. (2)). 

Again a convergence toward higher accuracy is seen in these results 
as the number of applied forces is increased. An important feature of the 
results is that the accuracy of the dependent DOF, calculated from Eg. 
(2-17), follows closely the same trends as the independent DOF. Also to be 
noted in the diagrams is the error levels for five shakers and for eight 
shakers. Although three shakers are exciting dependent DOF in the eight 
shaker case, the accuracy is either maintained or improved in almost every 
instance. Even if there happens to be a higher error indication for an 
increase in shakers, there is a definite lessening in the error spread. 

Two more series of tests were run. The first was a set of runs 
identical to the first test case with the exception that the forcing 
functions were swept-sines instead of random excitation. The results of 
this set were basically identical to the random excitation set. The final 
series also used the same set-up as the first test, except that no noise 
was introduced and the test range was from 0-1024 Hz. with 256 data points 
(a delta frequency of over 4 Hz.). The ninth mode corresponding to 25 Hz. 
was also deleted. All modes wen identified to extremely high accuracy 
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(the error was less than one part per billion) in spite of the fact that 
there were seven natural frequencies between two of the data points. 

Some preliminary testing of the algorithm utilizing the concept of 
pseudo-DOF was performed with mixed results. When a large percentage of 
the needed independent DOF were actually pseudo-DOF, there was a 
significant degradation in the accuracy of the results. However, if the 
pseudo-DOF were used sparingly to obtain the number of needed DOF, the 
algorithm performed on the same level as in the previously cited tests. In 
either case, the use of pseudo-DOF in spurious mode detection was 
satisfactory. An increase in the number of pseudo-DOF beyond the number of 
needed DOF had no ill effect on the algorithm's performance. In fact it 
helped facilitate the spurious mode identification. 
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Mode No. Frequency (Hz) Damping 


1 

5.000 

0.010 

2 

5.250 

0.010 

3 

5.500 

0.010 

4 

5.550 

0.010 

5 

5.900 

0,010 

6 

6.200 

0.010 

7 

6.600 

0.010 

8 

15.000 

0.010 

9 

25.000 

0.010 


DOF 1 
MODE 

2 

3 

4 

5 

6 

7 

8 

1 

.010 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

2 

1.000 

.010 

1.000 

1.000 

- 1.000 

- 1.000 

- 1.000 

- 1.000 

3 

1.000 

1.000 

-.010 

- 1.000 

- 1.000 

- 1.000 

1.000 

1.000 

4 

1.000 

1.000 

- 1.000 

-.010 

1.000 

1.000 

o 

o 

o 

1— f 

t 

- 1.000 

5 

1.000 

- 1.000 

- 1.000 

1.000 

.010 

- 1.000 

- 1.000 

1.000 

6 

1.000 

I 

1 — 1 

o 

o 

o 

- 1.000 

1.000 

- 1.000 

.010 

1.000 

-1.000 

7 

1.000 

- 1.000 

1.000 

- 1.000 

- 1.000 

1.000 

-.010 

1.000 

8 

1.000 

- 1.000 

1.000 

-1.000 

1.000 

- 1.000 

1.000 

-.010 

9 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 


y 

No. of Forces Force Locations I 

1 1 1 

4 1-4 

8 1-8 


Table 1. System Frequencies, Damping Factors, and Mode Shapes 

(defined at 8 DOF only); Forces and Locations ! 

'it 

H 
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Undamped Natural Frequency s .500630E+01 
Damping « .143159E-01 


Coordinate 

1 

2 

3 

4 

5 

6 

7 

8 


Coefficient 
• 308015E-02 
. 994932E+00 
.998240E+00 
.996091E+00 
.998748E+00 
• 998139E+00 
, 100000E+01 
.999306E+00 


-.277992E-01 

-.215889E-01 

-.101613E-01 

180890E-01 

-.492194E-02 

-.695993E-02 

0 

-.251417E-02 


Undamped Natural Frequency s .550004E+01 
Damping s . 1001.56E-01 


Coordinate 

1 

Coefficient 

•996047E+00 

-. 769268E-02 

2 

.996975E+00 

-.589179E-02 

3 

-.735022E-02 

.515352E-02 

4 

-.996900E+00 

.602848E-02 

5 

-.998680E+00 

.256406E-02 

6 

-.998711E+00 

.250478E-02 

7 

. 100000E+01 

0 

8 

•999892E+00 

-.209626E-03 


Undamped Natural Frequency - .554982E+01 
Damping = .999716E-02 


Coordinate 

1 

2 

3 

4 

5 

6 

7 

8 


Coefficient 
-.994959E+00 
-.992961E+00 
.998700E+00 
.954467E-02 
-.997289E+00 
-.997345E+00 
. 100000E+01 
.999725E+00 


-.443466E-03 
-.618205E-03 
. 114035E-03 
.389808E-04 
-.238917E-03 
-.234214E-03 
0 

.238590E-04 


No. of Forces = 1 


Noise Ratio » 0.0 


Table 2 
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Undamped Natural Frequency = .549B15E+01 
Damping * .922684E-02 


Coordinate 

1 

Coefficient 
. 100000E+01 

0 

2 

•961693E+00 

.154095E-01 

3 

~. 111589E+00 

-.243182E+00 

4 

~. 953094E+00 

189219E-01 

5 

-.866331E+00 

.201760E+00 

6 

-.874222E+00 

.225322E+00 

7 

.852086E+00 

-.247594E+00 

8 

.852671E+00 

-.202553E+00 


Undamped Natural Frequency - . 556574E+01 
Damping = .125684E-01 


Coordinate 

1 

2 

3 

4 

5 

6 

7 

8 


Coefficient 
.969798E+00 
. 100000E+01 
-.601954E+00 
-.306331E-01 
.663191E+00 
.655221E+00 
-.574999E+00 
-.641525E+00 


.718633E-02 

0 

.299185E+00 
. 115159E+00 
-. 152995E+00 
-. 18O12OE+O0 
.258952E+00 
. 178636E+00 


No. of Forces = 1 


Noise Ratio =0.1 


No mode corresponding to 5 Hz was identified 


Table 3 


} 
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Undamped Natural Frequency 3 .550027E+01 
Damping 3 .979364E-02 


Coordinate 

Coefficient 


1 

.897034E+00 

-.161891E-01 

2 

, 100000E+01 

0 

3 

-.263962E+00 

-.827695E-01 

4 

-.717776E+00 

.560193E-01 

5 

923468E-01 

.215160E+00 

6 

-. 173389E+00 

. 196507E+00 

7 

. 55 1836E+00 

-.126410E+00 

8 

.510949E+00 

- . 138541E+00 


Undamped Natural Frequency = .558794E+01 
Damping 3 .110378E-01 


Coordinate Coefficient 

1 -.721310E+00 

2 -.560135E+00 

3 .833390E+00 

4 -.769025E-01 

5 -.414973E+00 

6 -.508665E+00 

7 . lOOOOOE+Ol 

8 . 947 173E+00 


.649048E-01 
. 146199E-01 
-.662602E-02 
-.250925E-01 
122354E-01 
-.241696E-01 
0 

.470355E-02 


No. of Forces = 1 


Noise Ratio = 0.2 


No mode corresponding to 5 Hz was identified 


ssas 


Table 4 


Undamped Natural Frequency • .499677E+01 
Damping B .995318E-02 


Coordinate 

1 

1 ^ ’ 

Coefficient 

.140079E-01 

-.289103E-03 

2 

. 100000E+01 

0 

3 

.998304E+00 

. 122927E-03 

4 

.999471E+00 

• 385208E-04 

5 

♦997609E+00 

. 171764E-03 

6 

.997915E+00 

. 149431E-03 

7 

.996893E+00 

.22332^^-03 

8 

.997265E+00 

. 196417E-03 


Undamped Natural Frequency - .549998E+01 
Damping ■ . 100134E-0X 


Coordinate Coefficient 

1 . 1000Q0E+01 

2 .999560E+00 

3 112414E-01 

4 -.999586E+00 

5 -.998737E+00 

6 -.998722E+00 

7 . 9981 14E+00 

8 .998166E+00 


0 

. 154268E-02 
.434173E-02 
142288E-02 
-.438455E-02 
-.443570E-02 
• 657796E-02 
.639820E-02 


Undamped Natural Frequency = .554972E+01 
Damping = .999856E-02 


Coordinate 

1 

Coefficient 

-.992035E+00 

-.217490E-03 

2 

-.988883E+00 

-.304124E-03 

3 

.997943E+00 

. 5771 18E-04 

4 

•927480E-02 

.239258E-04 

5 

-.995724E+00 

-. 113431E-03 

6 

-.995813E+00 

-.110535E-03 

7 

. 100000E+01 

0 

8 

.999567E+00 

. 122911E-04 


No. of Forces - 4 


Noise Ratio - 0.0 
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Undamped Natural Frequency « .499588E+01 
Damping K . 101836E-01 


Coordinate 

Coefficient 


1 

.180826E-01 

.648566E-02 

2 

.989374E+00 

.422064E-03 

3 

.984854E+00 

.601076E-02 

4 

.lOQOOOE+Ol 

0 

5 

.975333E+00 

, 364940E-02 

6 

.987170E+00 

-,791028E-03 

7 

. 97755 XE-f 00 

.291086E-02 

8 

•975623E+00 

161876E-03 


Undamped Natural Frequency = .549835E+01 
Damping K .979900E-02 


Coordinate 

Coefficient 


1 

995591 E+00 

.353937E-01 

2 

-.963411E+00 

. 140834E-01 

3 

193424E-01 

- . 1 12736E-01 

4 

.987221E+00 

. 125764E-01 

5 

.952534E+0Q 

-.241061E-01 

6 

. 100000E+01 

0 

7 

995252E+00 

-,490844E-01 

8 

-.960618E+00 

- . 49 1432E-01 


Undamped Natural Frequency = .554825E+Q1 
Damping « .943557E-02 


Coordinate Coefficient 

1 . 993486E+00 

2 . 100000E+01 

3 -.973952E+00 

4 -.554699E-01 

5 .851456E+00 

6 . 903672E+00 

7 -.925752E+00 

8 -.905322E+00 


-.268232E-01 

0 

.261893E-02 

.508969E-02 

-.601461E-02 

198642E-02 

.291869E-01 

.464260E-01 


No. of Forces - 4 


Noise Ratio =0.1 


Table 6 
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I 


Undamped Natural Frequency - . 500344E+01 
Damping - .949072E-02 


Coordinate Coefficient 

1 -.431481E-02 

2 .981993E+00 

3 .924128E+00 

4 . 944983E+00 

5 . 97X504E+00 

6 , 100000E+01 

7 v 953668E+00 

8 .975622E+00 


**,383221E-02 
-.220043E-02 
, 172318E-02 
.305296E-01 
134243E-01 
0 

-.808437E-02 

-.367880E-02 


Undamped Natural Frequency a .548780E+01 
Damping 3 .959778E-02 


Coordinate 

Coefficient 


1 

.709776E+00 

. 119457E+00 

2 

.676447E+00 

. 170732E-01 

3 

, 129852E+00 

- . 441760E-01 

4 

-.915645E+00 

-.537032E-01 

5 

=.857429E+00 

«•. 136665E i r00 

6 

-.941759E+00 

145444E+00 

7 

. 100000E+01 

0 

8 

.955578E+00 

.678429E-01 


Undamped Natural Frequency = .554758E+01 
Damping 3 .104780E-01 


Coordinate 

1 

Coefficient 

.992678E+00 

.668940E-01 

2 

.100000E+01 

0 

3 

-.878481E+00 

-. 152797E+00 

4 

-.856727E-01 

. 457745E-01 

5 

.638683E+00 

. 176693E+00 

6 

.692493E+00 

. 154775E+00 

7 

- . 752810E+00 

-.214112E+00 

8 

-.713608E+00 

181675E+00 


4 


No. of Forces = 4 


Noise Ratio =0.2 


sirs#** 'W**r-«~c«rr,£r 


r; 

■i 


Table 7 
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Undamped Natural Frequency * .499655E+01 
Damping * ,985967E*02 


Coordinate 

1 

Coefficient 
, 142906E-01 

.864394E-03 

2 

.100000E+01 

0 

3 

•998184E+00 

-.367753E-03 

4 

.999433E+00 

-.U5383E-03 

5 

.997441E+00 

H.513542E-03 

6 

.997768E+00 

-.446875E-03 

7 

.996674E+00 

-.667799E-03 

8 

.997072E+00 

-.587373E-03 


Undamped Natural Frequency « .549997E+01 
Damping « ,999536E-02 


Coordinate 

1 

Coefficient 
. 100000E+01 

0 

2 

• 999341E+00 

.535210E-03 

3 

118561E-01 

.150597E-02 

4 

-.99P390E+00 

-.492155E-03 

K 

W 

^i998124E+Q0 

-.151908E-02 

6 

-.998102E+00 

153676E-02 

7 

.997187E+00 

.228076E-02 

8 

.997264E+00 

.221846E-02 


Undamped 

Natural Frequency 

« .554991E+01 

Damping - .999709E-02 

Coordinate 

1 

Coefficient 

-.997589E+00 

-.458045E-03 

2 

-.996632E+00 

-.638380E-03 

3 

.999379E+00 

. 116595E-03 

4 

.978273E-02 

.372708E-04 

5 

-.998702E+00 

-.249320E-03 

6 

-.998729E+00 

-.24478IE-03 

7 

. 100000E+01 

0 

8 

.999868E+00 

. 243785E-04 


No, of Forces * 8 Noise Ratio = 0.0 


Table 8 


Undamped Natural Frequency * . 499798E+01 
Damping « .965077E-02 


Coordinate 

1 

Coefficient 

.251566E-01 

-.728226E-02 

2 

•986205E+00 

-.423328E-02 

3 

.100000E+01 

0 

4 

.984109E+00 

513716E-02 

5 

,9734992+00 

.805894E-03 

6 

•971472E+00 

.548334E-02 

7 

•973073E+00 

.201598E-01 

8 

,9694752+00 

. 1419262-01 


Undamped Natural Frequency * .550583E+01 
Damping * •102792E-01 


Coordinate 

Coefficient 


1 

-.985845E+00 

- • X63209E-01 

2 

- . 994625E+00 

.324345E-01 

3 

.276250E-01 

-.977 544E-02 

4 

. 100000E+01 

0 

5 

.946089E+00 

.302677E-0] 

6 

. 952970E+00 

. 644X632-01 

7 

986397E+00 

. 187762E-02 

8 

-.959759E+00 

-.215829E-01 


Undamped Natural Frequency = .554904E+01 
Damping « .999550E-02 


Coordinate 

1 

2 

3 

4 

5 

6 

7 

8 


Coefficient 
.983620E+00 
. 100000E+01 
951970E+00 
-.645294E-01 
.875027E+00 
.874207E+00 
-.886737E+00 
-.890127E+00 


.708132E-02 

0 

. 111605E+00 
112379E+00 
191167E+00 
193172E+00 
.218127E+00 
. 182806E+00 


No. of Forces = 8 


Noise Ratio = 0. 


Table 9 
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Undamped Natural Frequency a .500069E+01 
Damping * .929434E-02 


Coordinate 

1 

2 

3 

4 

5 

6 

7 

8 


Coefficient 

.204826E-03 

♦995489E+00 

.988078E+00 

.995600E+00 

.995954E+00 

.980042E+00 

.985900E+00 

.lOQOOOE+Ol 


.548U1E-01 
. 149953E-01 
.501376E-01 
.579839E-01 
.971265E-Q2 
.484835E-02 
. 1 74225 E -02 
0 


Undamped Natural Frequency « .549181E+01 
Damping a .786745E-02 


Coordinate 

Coefficient 


1 

. 100000E+01 

0 

2 

.978124E+00 

.865004E-01 

3 

-.562677E-02 

.274473E-01 

4 

-.949020E+00 

-.470651E-01 

5 

-.928190E+00 

-. 187158E+00 

6 

-. 959287E+00 

-.153123E+00 

7 

.969323E+00 

.904357E-01 

8 

.982434E+00 

. 148072E+00 


Undamped Natural Frequency - .555115E+01 
Damping ~ .100090E-01 


Coordinate 

1 

2 

3 

4 

5 

6 

7 

8 


Coefficient 
-.835000E+00 
-.823659E+00 
.921431E+00 
-.618963E-01 
-.966994E+00 
-.963982E+00 
, 100000E+01 
.969904E+00 


-. 188699E+00 
-. 144876E+G0 
. 105953E+00 
.955734E-01 
-.642726E-01 
-.450693E-01 
0 

.349342E-01 


No. of Forces = 8 


Noise Ratio = 0.2 


Table 10 
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Mode No. 

Frequency (Hz) 

Damping 

1 

5.000 

0.010 

2 

5.250 

0.010 

3 

5.500 

0.010 

4 

5.550 

0.010 

5 

15.000 

0.010 

6 

25.000 

0.010 


DOF 1 
MODE 

2 

3 

4 

5 

6 

7 

8 

1 

.010 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

2 

1.000 

.010 

1.000 

1.000 

- 1.000 

- 1.000 

- 1.000 

- 1.000 

3 

1.000 

1.000 

-.010 

- 1.000 

o 

o 

o 

1 — 1 
t 

- 1.000 

1.000 

1.000 

4 

1.000 

1.000 

- 1.000 

-.010 

1.000 

1.000 

- 1.000 

- 1.000 

5 

1.000 

- 1.000 

- 1.000 

1.000 

.010 

- 1.000 

- 1.000 

1.000 

6 

1.000 

- 1.000 

- 1.000 

1.000 

- 1.000 

,010 

1.000 

- 1.000 


No. of Forces 
1 
5 
8 


Force Locations 
1 

1-5 

1-8 


Table 11. System Frequencies, Damping Factors, and Mode Shapes 
(defined at 8 DOF only); Forces and Locations 












Chapter 4 


4.1 Conclusions 

In this thesis, the theory of a method used to determine an estima- 
tion of the modal parameters for a viscously-damped, linear structure has 
been developed. The MSMA method involves the application and measurement 
of several non-coherent forces. Verification of the method was 
demonstrated analytically through an extensive computer simulation study. 

The computer simulation showed the MSMA method to be stable. With 
increasing amounts of noise the results did not diverge dramatically. 
Also, as the number of applied forces was increased, the results converged 
to the system values. 

Several advantages of the MSMA method were demonstrated by the 
computer analysis. The first is the ability to locate a mode which has a 
node line in the vicinity of a shaker by applying several forces. The 
second is the small amount of user judgement and expertise involved. 
Basically the only user judgement needed is a knowledge of the approximate 
number of modes in the analysis range, and that is usually supplied by a 
previous finite element simulation or can be estimated by examination of a 
preliminary FRF. The third advantage is the relatively fast analysis time. 
With the small amount of user interaction the analysis time is basically 
dependent upon computation speed and total amount of data processed. 
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4.2 Recommendations for Further Work 


There are still many features of the MSMA method that need to be 
examined. 

a. A modal confidence factor (MCF) should be developed to system- 
atically evaluate the validity of acquired modes. The MCF should 
be based on a comparison of the eigenvectors corresponding to the 
pseudo-DOF and the eigenvectors of their companion true-DOF. 

b. A preliminary mode count procedure should be established. 
Instead of an initial count based on a finite element model, there 
should be a method of determining an estimate from the structure 
itself. One method may be to administer SPR tests at a couple of 
points to obtain a rough estimate, and then increase that by some 
margin. 

c. To make this method more efficient and more capable of being 
put into small computers, numerical methods to reduce the size of 
the data block and required working space should be developed. 

d. Since the [C .] and [K..] matrices are system constants, it 
might prove advantageous to use cross-spectrums of the forces and 
accelerations to reduce the statistical variance. One form of this 
technique would modify Eq. (2-9) to obtain 


-CG(XX) G ( XX ) G(X>)] = [tC 1 ][K 1 .]CD. d ]] 


G ( XX ) 

G(XX) 

G(XF) 

G(XX) 

G(XX) 

G{ X F) 

g(fJ) 

G(FX) 

G(FF) 

A and B. 

This 

equation 

matrices 

in one of two 
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I* 


fashions. The first would be to average the cross-spectrum 
matrices for each frequency point into one matrix before carrying 
out the mathematical operations and subsequent solution for the 
modal matrices. The second would be to obtain an estimate for the 
modal matrices for each frequency point and then average those 
estimates together to obtain the overall estimate. 

e. Though the equations from Chapter 2 are applicable to any M, 
C, and K matrices, the derivation of the equations was for a system 
with proportional damping. A formal derivation for non- 
proportional damping should be developed. 

f. Alternative methods for the development of pseudo-DOF should 
be explored and developed. A method which uses data previously 
acquired and manipulates the data numerically, through a phasing 
scheme or some other approach, would be highly attractive in that 
it would eliminate additional testing of the structure to produce 
the pseudo-DOF. 


APPENDIX 


Computer Program Listing 


PROGRAM RCMB2(TTY , INPUT .OUTPUT ,TAPE2=TTY ,TAPE5=INPUT , TAPE6=OUTPUT , 
1TAPE1) 


IN THE PROGRAM ALL INTEGERS ARE INPUT IN 15 FORMAT, AND ALL 
REAL NUMBERS IN F10 FORMAT. 

THE USER IS FIRST ASKED HOW MANY D.O.F. THE SYSTEM HAS( I.F.. 

HOW MANY RESPONSE LOCATIONS). 

NEXT THE USER IS ASKED HOW MANY MODES ARE TO BE IN THE SYSTEM 
(THE NUMBER OF MODES DO NOT HAVE TO EQUAL THE NUMBER OF D.O.F.). 

THE FREQUENCIES AND DAMPING RATIOS ARE THEN INPUT. 

THE MODE SHAPES ARE INPUT AFTER THAT. 

THE NEXT INPUT IS THE INITIAL FREQUENCY, AND THE FINAL FREQUENCY 
GIVING THE FREQUENCY RANGE TO BE TESTED IN 256 INCREMENTS. 

THE PROGRAM WILL OUTPUT THE FREQUENCIES CORRESPONDING TO THE 
MAXIMUM AND MINIMUM AMPLITUDES IN THE IMAGINARY RESPONSE. 

THE USER HAS THE CHOICE OF CHOOSING A DIFFERENT RESPONSE POINT 
(THE REFERENCE POINT IS FIXED AT 1) OR JUST USING THE GIVEN 
VALUES FOR ESTIMATES. 

THE USER IS THEN ASKED HOW MANY MODES ARE TO BE EXAMINED( THIS 
CANNOT EXCEED THE NUMBER OF D.O.F.), AND IS THEN REQUESTED TO 
INPUT THE CORRESPONDING NUMBER OF RESPONSE LOCATIONS, AND THE 
EXCITATION LOCATION. 

THE PROGRAM RESPONDS .WITH THE EIGENVALUES AND EIGEN VECTORS. 

THIS PROGRAM USES IMSLIBF ROUTINE 1 EIGRF' 

AND 'GGUBFS' 


C0MM0N/FREQ/AM(20) ,W(20) ,D(20) ,SM(20,20) .IFOR(IO) ,NF 
COMMON/DATA/AI(256) , AR( 256 ) ,WW(256) ,F(10,512) ,UU(256) 
COMMON/PA/S(30.1O) ,SMBKI(30,30) ,SMBK(10,30) .ASIN(IO.IO) 
C0MM0N/AR/XX(30,30) .SINT(IO.IO) .SMINV(IOO) .SMINVT(IOO) 
DIMENSION Y(10,512) ,X(30,512) ,YX(30,30) ,XIN(30,30) 
DIMENSION T(10,30) .IRESRE(IO) 

DIMENSION IDL(10),YD(10,512) 

DIMENSION ITYPE(20),IPHI(20),IRESL0C(20) 
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DIMENSION ZZZ(IOO) ,LB(29) ,MB(29) ,SY(400) 

DIMENSION IJK(20) » PDPI( 20 , XO) 

DIMENSION PMAG( 20 , 20) ,THETA(20,20) ,WK( 1000), XMEN( 900) 
COMPLEX AW(20) ,PHI(20 ,20) ,ZZ(20 ,20) 

DOUBLE PRECISION DSEED 

* SETTING CONSTANTS 
DSEED=1 .DO 
TW0RI=8.*ATAN(1.) 

* INPUTTING DATA 
WRITE( 2 , 80) 

PRINT 80 

80 FORMAT (IX, 'ARE YOUR INPUT FREQUENCIES AND SHAPES ON A DATA 
1/, IX, 'FILE? [Y=1:N=Q]' ) 

PRINT 27 
WRITE(2,27) 

READ(2, 400)1 
PRINT 301,1 
IF(I.NQ.l) GOTO 600 
PRINT 1 

READ(5,400)NDOF 
PRINT 301.NDOF 
PRINT 2 

READ(5 .4QQ)NFR 
PRINT 301, NFR 
PRINT 3 

DO 200 1=1, NFR 
PRINT 4,1 

READ(5,401) ,W(I) ,D(I) 

PRINT 300,W(I) ,D(I) 

W(I)*W(I) w TWOPI 

200 CONTINUE 
PRINT 5 
INT=ND0F/7+l 
IEXTRA=ND0F-(INT“1)*7 
DO 201 1=1, NFR 
PRINT 6,1 

DO 201 JJ=1 ,INT 
IEND=7 

IF(JJ.EQ.INT) IEND=IEXTRA 

IF(IEND.EQ.O) GOTO 201 

READ(5 ,401)(SM((JJ-1)*7+J ,1) ,J=1 ,IEND) 

201 CONTINUE 

DO 202 1=1, NFR 

202 PRINT 300,(SM(J,I) ,J=1,ND0F) 

GOTO 601 

600 CONTINUE 
WRITE( 2 , 14) 

PRINT 14 

14 F0RMAT(1X,'Y0U MAY ENTER UP TO 20 MODES AND D.O.F. AND', 
1/, IX, ' ANALIZE UP TO 10.',//) 


PRINT 1 
WRITE(2,1) 

PRINT 27 
WRITER 27 ) 

27 FORMAT(ix, [ 15 FORMAT ]') 

28 F0RMAT(1X, 1 [ F10 FORMAT ]’) 

1 FORMAT( IX, 'ENTER NUMBER OF DEGREES OF FREEDOM') 
READ(2,40Q) NDOF 

400 FORMAT (15) 

PRINT 301, NDOF 

300 F0RMAT(1X, ' *' ,19(1X,F6.3)) 

301 FORMAT ( IX , ,2X,4(2X,I5)) 

PRINT 2 

WRITE(2 ,2) 

PRINT 27 
WRITE(2,27) 

2 FORMAT ( IX , 1 ENTER NUMBER OF MODES ' ) 

READ(2,40Q) NFR 

PRINT 301 , NFR 
PRINT 3 
WRITE(2,3) 

PRINT 28 
WRITE(2,28) 

3 FORMAT( IX, 'ENTER FREQUENCIES [IN HZ,], DAMPING RATIOS') 
DO 100 1=1, NFR 

PRINT 4,1 
WRITE(2,4)I 

4 F0RMAT(1X,' FREQUENCY NO.', 13) 

READ(2,401) W(I),D(I) 

401 FORMAT (7(E10 .4)) 

PRINT 300 , W( I ) » D( I ) 

W(I)=TWOPI*W(I) 

100 CONTINUE 
PRINT 5 
WRITE(2,5) 

PRINT 28 
WRITE(2,28) 

5 FORMAT( IX, ' ENTER MODE SHAPES ' ) 

INT=ND0F/7+l 

IEXTRA=ND0F-(INT~1)*7 

DO 101 1=1, NFR 
PRINT 6,1 
WRITE( 2 ,6)1 

6 FORMAT( IX,' ENTER MODE NO.', 13,' SEVEN D.O.F. PER LINE') 
DO 101 JJ=1 , I NT 

IEND=7 

IF(JJ.EQ.INT) IEND=IEXTRA 

IF(IEND.EQ.O) GOTO 101 

READC 2 , 401 ) ( SM( ( J J-l )*7+J , I ) , J=1 , 1 END) 

101 CONTINUE 
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DO 102 I-l.NFR 

102 PRINT 300,(SM(J,I) ,J=1 ,NDOF) 

601 CONTINUE 

DO 569 1=1, NFR 
AM( I )— 0 . 0 
DO 569 J=1 , NDOF 

569 AM(I)=AM(I)+SM(J,I)*SM(J,I) 

888 CONTINUE 

SETTING PARAMETERS FOR RUN 

WRITE(2,7) 

PRINT 7 
PRINT 28 
WRITE(2,28) 

7 FORMAT( IX, 'ENTER INITIAL FREQUENCY RANGE IN HZ. ; INIT. FREQ. / ' 
1* FINAL FREQ,') 

READ(2,401) STFR.FIFR 
PRINT 330, STFR.FIFR 
330 FORMAT ( IX , 1 *' ,2(1X,F12.5)) 

STFR=TWOPI*STFR 
FIFR=TWOPI*FIFR 
FIST=FI FR-STFR 
WRITE(2,29) 

PRINT 29 
WRITE(2,27) 

PRINT 27 

READ(2,400) ITWO 
PRINT 301, ITWO 

29 F0RMAT(1X,' INPUT THE NUMBER OF FREQUENCY POINTS TO BE ANALIZED' , 
1/,1X,'[UP TO 00256]') 

IDTWO=2*ITWO 
ITT=ITWO-l 
TW0M1= FLOAT (ITT) 

WRITE(2,32) 

PRINT 32 
WRITE( 2 , 28) 

PRINT 28 

READ(2,401) ANOISE 
PRINT 300, ANOISE 

32 F0RMAT(1X,' INPUT THE PERCENT OF NOISE IN THE FUNCTIONS', 

1/,1X,'[ 100% =!.]') 

GENERATING A CLEAN TRANSFER FUNCTION REFERENCED TO DOF 1. 

PRINTS OUT ALL LOCAL MAXIMUMS AND MINIMUMS IN THE IMAGINARY 
PUNE. 

519 WRITE( 2,16) 

PRINT 16 
°RINT 27 


f. 


WRITE(2 27) 

16 FORMAT(ix, ' ENTER THE RESPONSE POINT FOR THE INITIAL RUN') 
READ( 2 , 400) JRESP 
PRINT 301, JRESP 


INITIALIZE PARAMETERS TO KEEP FROM GETTING A WARNING 


ANSOLD=0. 

DIFOLD=0, 

WWQLD=Q, 

DO 106 J=1 , I TWO 
TT"FLOAT(J) 

WW( J )-STFR+( TT- 1 . ) * FI ST/TWOM1 
CALL TRAN(WW(J) ,1,ANS,NFR, JRESP) 

IF(J.EQ.l) GOTO 123 
DIF=ANS-ANSOLD 
IF(J.EQ.2) GOTO 124 
POS=DI F^DI FOLD 
IF(POS.GE.O) GOTO 124 
WRITE(2,12)WW0LD,ANS0LD 
PRINT 12,WW0LD,ANS0LD 
124 DIFOLD=DIF 
123 ANSOLD=ANS 

WWO LD=WW( J ) /TWOP I 

12 FORMAT ( IX , l W= l ,E12,6, ' AMP=',E12.6) 

106 CONTINUE 
WRITE(2, 17) 

PRINT 17 
PRINT 27 
WRITE(2,27) 

17 F0RMAT(1X, 1 IS ANOTHER RESPONSE POINT REQUIRED? [Y*1/N=0]') 
READ( 2,400) IANS 
PRINT 301, IANS 
IF( IANS.EQ. 1) GOTO 519 
WRITE(2,8) 

8 F0RMAT(1X> 5 ENTER NUMBER OF MODES TO BE EXAMINED') 

PRINT 8 
PRINT 27 
WRITE(2,27) 

READ(2,400) NME 
PRINT 301, NME 
WRITE(2, 10)NME 
PRINT 10, NME 
PRINT 27 
WRITE(2,27) 

10 FORMATC IX,' ENTER THE 1 ,13, 

1' DIFFERENT INDEPENDENT RESPONSE LOCATIONS;') 

1/, IX, 'ONE PER LINE 1 ) 


54 
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DO 103 1=1, NME 
READ(2,400) IRESLOC(I) 

PRINT 301 ,IRESLOC(I) 

103 CONTINUE 
WRITE(2 , 26) 

PRINT 26 
WRITE( 2,27) 

PRINT 27 

26 FORMAT( IX, 'ENTER THE NUMBER OF FORCED INDEPENDENT RESPONSE 1 , 
1/ , IX LOCATIONS') 

READ(2 ,400)NF 
PRINT 301 „NF 
DO 509 1=1, NF 
WRITE( 2,40) 

PRINT 40 

READ(2 ,400)IFOR( I) 

PRINT 301 ,IFOR(I) 

509 CONTINUE 

40 FORMAT( IX,' INPUT FORCE LOCATION') 

WRITE( 2,21) 

PRINT 21 
PRINT 27 
WRITE( 2,27) 

21 FORMATC IX , ' ENTER THE NUMBER OF DEPENDENT RESPONSE LOCATIONS') 
READ(2 ,400)NDL 

PRINT 301 , NDL 
IF(NDL.EQ.O) GOTO 505 
WRXTE( 2 , 22 )NDL 
PRINT 22, NDL 
PRINT 27 
WRITE(2,27) 

22 FORMAT( ix, ' ENTER THE* ,13, ' DEPENDENT RESPONSE LOCATIONS;', 

1/, IX, 'ONE PER LINE') 

DO 522 1=1, NDL 
READ(2 ,400)IDL(I) 

522 PRINT 301,IDL(I) 

505 CONTINUE 
603 CONTINUE 
WRITE(2,602) 

PRINT 602 
WRITE(2,27) 

PRINT 27 

READ( 2 , 400) IQUEST 
PRINT 301, IQUEST 

602 FORMAT( IX, 'DO YOU WANT SWEPT SINE [1]', 

1/, 'OR RANDOM [2] FORCING?') 

IF(IQUEST.NE. 1. AND, IQUEST. NE. 2) GOTO 603 
IF( IQUEST. EQ. 2) GOTO 605 
WRITE( 2,606) , NF 
PRINT 606, NF 


DO 604 1=1, NF 
WRITE(2,28) 

READ(2,401)LB(I) 

PRINT 300,LB(I) 

604 CONTINUE 

606 FORMAT ( XX , 1 ENTER THE* ,12,' PHASE SHIFTS') 
WR1TE(2,608) 

PRINT 608 
WRITE(2,27) 

PRINT 27 

READ(2,400) IAVE 
PRINT 301, IAVE 

608 FORMAT ( IX » 1 HOW MANY AVERAGES DO YOU WANT?') 

DO 576 KK=1,IAVE 
DO 576 1=1 ,NDOF 
DO 577 J=1,NF 
K=IFOR(J) 

577 IF(K.EQ.I) GOTO 578 
DO 579 J=l, IDTWO 

579 F( I ,J)=0.0 
GOTO 580 

578 CONTINUE 
PHII S LB( J) 

CALL ASINE( FIST, IDTWO, STFR , FI FR , PH I I , TW0M1 , XMEN ) 
DO 582 11=1, IDTWO 
582 Z(II)=CMPLX(XMEN(II) ,0.0) 

CALL FFT( IDTWO, Z,TWOPI) 

DO 581 J=1 , ITWO 
UJJ=J+(KK~1)*ITW0 
UJ=JJJ+IAVE*ITWO 
F(K,JJJ)=REAL(Z(J)) 

F(K,JJ)=AIMAG(Z(J)) 

581 CONTINUE 

580 CONTINUE 
576 CONTINUE 

GOTO 607 

605 CONTINUE 
WRITE(2,608) 

PRINT 608 
WRITE(2,27) 

PRINT 27 

READ(2,400) IAVE 
PRINT 301, IAVE 
DO 568 KK=1,IAVE 
DO 568 I=1,ND0F 
DO 570 J=1,NF 
K=IFOR(J) 

570 IF(K.EQ.I) GOTO 571 
DO 572 J=l, IDTWO 
572 F(I ,J)=0.0 
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GOTO 573 
571 CONTINUE 

DO 574 J=l,ITWO 
JJJ=J+(KK-l)*ITWO 
JU=JJJ+IAVE*ITWO 
ARAN=GGUBFS(DSEED) 

ARAN=ARAN*TWOPI 

F(K,JJJ)=COS(ARAN) 

F(K,JJ)=SIN(ARAN) 

574 CONTINUE 
573 CONTINUE 
568 CONTINUE 
607 CONTINUE 

* 

* Y IS AN (NME)*(512) MATRIX CONTAINING THE ACCELERATION RESPONSES 

* THE REAL PARTS ARE IN TNE FIRST 256 COLUMNS AND THE IMAGINARY ARE 

* * IN THE SECOND 256 COLUMNS. 

* 

* X IS AN (2*NME+1)*(512) MATRIX CONTAINING THE VELOCITY RESPONSES 

* IN THE FIRST NME ROWS ( IN THE SAME FASHION AS THE ACCELERATION 

* RESPONSES), THE DISPLACEMENT RESPONSES IN THE SECOND NME ROWS. 

* THE FINAL ROW CONTAINS 1. IN THE FIRST 256 COLUMNS AND 0. IN 

* THE FINAL 256. 

* 

* 

DO 500 K— 1 , IAVE 
DO 500 1=1, NME 
IN=IRESLOC(I) 

CALL FRF( ITWO » TW0M1 , STFR , FI FR , NFR , IN , ANOISE , DSEED , K , IAVE) 

DO 500 J=1 ,ITWO 

JJJ=J+(K~1)*ITW0 

JJ=JJJ+IAVE*m/0 

II=I+NME 

III=II+NME 

X(II ,JJ)=AI(J) 

X(II,JJJ)=AR(J) 

X(I ,JJJ)=-WW(J)*X(II ,JJ> 

X(I ,JJ)=V/W(J) A X(II ,JJJ) 

Y(I ,JJJ)=~WW(J)*WW(J)*X(II ,JJJ) 

Y(I ,JJ)=-WW(J) yf WW(J)*X(II ,JJ) 

500 CONTINUE 

IDTAVE=IDTWO*IAVE 
DO 575 1=1, NF 
II-I+NME+NME 
K=IFOR(I) 

DO 575 J=1 ; IDTAVE 
X(II ,J)=F(K.J) 

575 CONTINUE 
IF(NDL.EQ.O) GOTO 530 
DO 523 K=1 ,IAVE 
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DO 523 I=1,NDL 
IN=IDL( I ) 

CALL FRF(ITWO,TWOMl,STFR,FIFR,NFR,IN,ANOISE,DSEED,K,IAVE) 

DO 523 J=l, ITWO 

JJJ=J+(K-l)*ITWO 

JJ=JJJ+IAVE*ITWO 

YD(I,JJJ)=“WW(J)*WW(J) A AR(J) 

YD( I , JJ)=-WW(U)*WW( J)*AI( J) 

523 CONTINUE 
530 CONTINUE 
N3~2*NME+NF 
N2-2*NME 
ND=NME-NF 
ND3=3*ND 
N F3==3 A NF 


* XIN=(X)*(TRANSPOSE(X)) 

* XIN IS PUT INTO AN ONE DIMENSIONAL ARRAY XMEN SO THAT IT CAN BE 

* INVERTED EASILY, THE INVERTED XMEN IS THEN PUT BACK INTO KIN 


DO 501 1=1, N3 
DO 501 J=1 ,N3 
XIN(I,J)=0. 

DO 501 K=1,IDTAVE 
XIN(I,J)=XIN(I,J)+X(I,K)*X(J,K) 
501 CONTINUE 
13 FORMAT( IX, 17(E8. 2)) 

ID=0 

DO 510 J=1,N3 
DO 510 1=1, N3 
ID=ID+1 

XMEN(ID)=XIN( I ,J) 

510 CONTINUE 

CALL MINV( XMEN, N3, DETER, LB, MB) 
ID=0 

DO 511 J=1 ,N3 
DO 511 1=1, N3 
ID=ID+1 

XIN(I ,J)=XMEN( ID) 

511 CONTINUE 


* YX=(Y)*(TRANSPOSE(X)) 

* T=(YX)*(XIN) ; XIN BEING THE INVERTED FORM 

* 

* 


DO 502 1=1, NME 
DO 502 J=1,N3 
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YX(I,J)=0, 

DO 502 K=1 , IDTAVE 
m YX(I ,J)=YX(I ,J)+Y(I,K)*X(J,K) 

DO 503 1=1, NME 
DO 503 J=1 ,N3 
T(I,J)=0. 

DO 503 K=1 ,N3 

503 T(I ,J)=T(I ,J)-YX(I ,K)*XIN(K,J) 

a 

* THE T MATRIX IS COMPOSED OF THREE SUBMATRICES; 

* THE FIRST IS B(NME,NME), THE NEXT IS K(NME,NME), THE 

* LAST IS GAMMA(NME,1). 

A 

* SY IS A MATRIX COMPOSED OF B,K,AN IDENTITY MATRIX, AND A 

* NULL MATRIX, AS FOLLOWS: 

A 

* SY = [ B K ] 

* [-1 0 ] 

A 

* THE NEGATIVE OF SY IS PUT INTO AN EIGEN SOLVER AND THE VALUES 

* RETURNED ARE OUTPUT. 

A 

* THE EIGENVALUES AND VECTORS ARE ALSO GENTERATED USING ONLY THE 

* K MATRIX. 

A 

A 

IT=-NME 

DO 504 J=1,2*NME 
IT=IT+NME 
DO 504 1=1, NME 
IT=IT+1 

504 $Y(IT)=T(I,J) 

IT=0 

N2=2*NME 

DO 506 1=1, N2 

IT=IT+NME 

DO 506 J=1 ,NME 

IT=IT+1 

SY(IT)=0. 

IF(I.EQ.J) SY(IT)=-1. 

506 CONTINUE 
514 CONTINUE 

DO 518 1=1 ,N2*N2 
SY(I)=-SY(I) 

518 CONTINUE 

CALL EIGRF(SY,N2,N2,2,AW,Z,N2,WK,IER) 

A 

* SCALING EIGENVECTORS 

it 

1 R— 0 
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DO 531 J=1,N2 
DO 531 1=1, N2 
IR=IR+1 
ZZ(I ,J)=Z(IR) 

531 CONTINUE 

DO 545 J=1 ,N2 
ZMAG0=0 . 0 
DO 545 11=1, NME 
I=II+NME 

ZMAG=CABS(ZZ(I,J)) 

PMAG(I ,J)=ZMAG 
IF( II .GT.NME) GOTO 42 
I F(ABS( PMAG( I ,J)).GT. ZMAGO ) I JK( J )=I 
I F(ABS( PMAG( I ,J)).GT. ZMAGO) ZMAGO=ABS(PMAG( I , J ) ) 
42 CONTINUE 

RE=REAL(ZZ( I ,J)) 

AII=AIMAG(ZZ( I , J ) ) 

THETA( I ,J)=ATAN2(AII,RE) 

545 CONTINUE 

DO 546 J=1 ,N2 
DO 546 11=1, NME 
I=II+NME 
K=IJK(J) 

THET=THETA( I , J)-THETA(K.J) 

RE=( PMAG( I , J ) *COS( THET ) )/ABS( PMAG( K , J ) ) 
AII=(PMAG(I ,J)*SIN(THET))/ABS(PMAG(K,J)) 

ZZ(I ,J)=CMPLX(RE,AII) 

546 CONTINUE 


* TO FIND THE MODE SHAPE COEFFICIENTS FOR THE DEPENDENT D.O.F. 

* USING THE EQUATION: 

* 

* PHI(DEP)*INV[PHI(IND)]=[Y*TRANS(X)]*INV[X*TRANS(X)] 

* WHERE [Y]= THE ACCELERATION RESPONSES FOR THE DEPENDENT D.O.F. 

* [X]= THE ACCELERATION RESPONSES FOR THE INDEP. D.O.F. 

* 

* PDPI= YX* YY WHERE YY HAS BEEN INVERTED 


IF(NDL.EQ.O) GOTO 529 
DO 524 1=1, NME 
DO 524 0=1, NME 
XIN(I ,J)=0.0 
DO 524 K=1 , IDTAVE 

524 XIN(I ,J)=XIN(I,J)+Y(I ,K)*Y(J,K) 
DO 525 1=1, NDL 
DO 525 J=1 ,NME 
YX( I , J )=0 . 0 
DO 525 K=l, IDTAVE 
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525 YX( I ,J)=YX< I , J )+YD( I , K)*Y(J , K) 

IT=0 

DO 526 J=1,NME 
DO 526 1=1, NME 
IT=IT+1 

526 XMEN(IT)=XIN(I,J) 

CALL MINV(XMEN , NME , DETER , LB ,MB) 

IT=0 

DO 527 J=1 , NME 
DO 527 1=1, NME 
IT=IT+1 

527 XIN(I ,J)=XMEN(IT) 

DO 528 1=1, NDL 

DO 528 J=1 ,NME 
PDPI(I,J)=0.0 
DO 528 K=1 ,NME 

528 PDPI(I,J)=PDPI(I,J)+YX(I,K)*XIN(K,J) 
ND2=2*NDL 

DO 532 1=1, NDL 
DO 532 J=1 ,N2 
PHI(I ,J)=0.0 
DO 532 K=1,NME 
KK=K+NME 

532 PHI(I ,J)=PHI(I ,J)+PDPI(I ,K)*ZZ(KK,J) 

529 CONTINUE 


* PRINT OUT THE CONVERGENCE CRITERIA FOR THE EIGEN SOLUTION 

* ( LESS THAN ONE INDICATES GOOD CONVERGENCE) 

WRITE(2 ,11)WK(1) 

PRINT 11 ,WK(1) 

* SORTING THE DEPENDENT AND INDEPENDENT D.O.F. INTO 

* CONSECTIVE ORDER 


ITOT=NDL+NME 
DO 533 1=1 ,ITOT 
DO 534 J=1 ,NME 
IF(IRESLOC(J).NQ.I) GOTO 534 
IPHI(I)=J 
ITYPE(I)=1 
GOTO 533 
534 CONTINUE 

IF(NDL.EQ.O) GOTO 533 
DO 541 J=1 ,NDL 
IF(IDL(J) .NQ. I) GOTO 541 
IPHI(I)=J 
ITYPE(I)=2 
541 CONTINUE 
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533 CONTINUE 


* PRINTING OUT THE EIGEN SOLUTIONS 


DO 547 K=1 ,N2 
AW( K )=AW( K )/TWOP I 
SIG=-REAL(AW(K)) 

DEL=AIMAG(AW(K)) 

SQ=SQRT(SIG*SIG+DEL*DEL) 

DMP=SIG/SQ 
PRINT 24,K,SQ,DMP 

24 FORMAT(/// f IX , 'MODE NO. M3, 2X, 'NATURAL UNDAMPED FREQUENCY= 
1E12. 6, /,33X, 'DAMPING = ' ,E12.6, 

1// , IX , ' COORDINATE' ,9X , ' COEFFICIENT' ,/) 

DO 535 1=1 , ITOT 
J=IPHI(I) 

JJ=J+NME 

I F( ITYPE( I ) . EQ . 2 ) GOTO 536 
PRINT 23,1 ,ZZ(JJ,K) 

GOTO 535 

536 PRINT 23,1 ,PHI(J,K) 

535 CONTINUE 

IF(NMEID.EQ.G) GOTO 43 
DO 44 1=1 , NMEID 
II=I+NME+NME 

44 PRINT 23,1 RESRE( I ) ,ZZ( I I , K) 

43 CONTINUE 
547 CONTINUE 

23 FORMAT(4X,I2,12X,E12.6,2X,E12.6) 


* SETTING THE EIGEN PROBLEM FOR THE STIFFNESS MATRIX 

* XMEN(NME,NME) = T(NME,[NME TO 2*NME]) 

IT=0 


DO 515 J=1 ,NME 
DO 515 1=1, NME 
JJ=J+NME 
IT=IT+1 

XMEN( IT)=T ( I ,JJ) 

515 CONTINUE 
PRINT 31 

31 F0RMAT(//,1X,'K MATRIX') 

DO 516 1=1, NME 

PRINT 18,(XMEN(I+(J-1)*NME) ,J=1 ,NME) 
18 FORMAT(3X ,10E12 .6) 

516 CONTINUE 




CALL EIGRF( XMEN , NME , NME , 2 , AW , Z , NME ,WK , IHR) 

IR=0 

DO 537 J=1 ,NME 
DO 537 1=1, NME 
IR=IR+1 
ZZ(I ,J)=Z(IR) 

537 CONTINUE 

DO 555 J=1,NME 

ZMA60=0 . 0 

DO 555 11=1 , NME 

PMAG( I I , J )=REAL( ZZ( II , J ) ) 

IF(II.GT.NME) GOTO 45 
IF(ABS(PMAG(II,J)).GT.ZMAGO) IJK(J)=II 
I F( ABS(PMAG( I I , J ) ) . GT . ZMAGO ) ZMAGO=ABS(PMAG( I I , J ) ) 
45 CONTINUE 

555 CONTINUE 

DO 556 0=1 , NME 
DO 556 11=1, NME 
K=IJK(J) 

ZZ(II,J)=ZZ(II,J)/ABS(PMAG(K,J)) 

556 CONTINUE 


* FINDING THE MODE SHAPE COEFFICIENTS FOR THE DEP. D.O.F. 


IF(NDL.EQ.O) GOTO 544 
DO 538 1=1, NDL 
DO 538 J=1 ,NME 
PHI ( I , J )=0 . 0 
DO 538 K=1,NME 

PHI(I ,J)=PHI(I , J)+PDPI(I ,K)*ZZ(K,J) 
538 CONTINUE 
544 CONTINUE 


A COMPUTING THE DAMPING FACTORS 

* DIAG[INV[PHI( IND)]*B(MATRIX)*PHI(IND)] =[2*DAMP( I)*AW(I)] 

* 


IT=0 

DO 508 J=1 ,NME 
DO 508 1=1, NME 
IT=IT+1 
ZZZ( IT)=Z( IT) 
PDPI ( I ,0 )=Z( IT) 
DO 542 1=1, NME 
DO 542 J=1 ,NME 
XIN(I,J)=0.0 
DO 542 K=1,NME 
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542 XIN(I ,J)=XIN(I ,J)+T(I ,K)WI(K ( J) 

CALL MINV( ZZZ , NME , DETER , LB , MB) 

IT=0 

DO 543 J=1,NME 
DO 543 1=1, NME 
IT=IT+1 

543 PDPI(I ,J)=ZZZ(IT) 

DO 520 1=1, NME 

DO 520 J=1,NME 
THETA(I ,J)=0.0 
DO 520 K=1 , NME 

520 THETA(I ,J)=THETA(I ,J)+PDPI(I ,K)*XIN(K,J) 


* CHANGING EIGENVALUES TO HERTZ AND OUTPUTTING SOLUTIONS 


DO 548 K=1 , NME 
AW(K)=AW(K)/(TWOPI*TWOPI) 

RAW=AW(K) 

ICOUNT=0 

IF( RAW. LT.O) IC0UNT=1 
RAW=ABS( RAW) 

RAW=SQRT(RAW) 

DAMP=THETA(K,K)/(2.*TW0PI*RAW) 

PRINT 20, K, RAW, DAMP 
IF(ICOUNT.EQ.l) PRINT 25 
25 FORMAT( IX, 'OMEGA SQUARED WAS NEGATIVE',//) 

20 F0RMAT(1X, 'MODE NO. ',12,' NATURAL FREQ.(HZ)= ',E12.6, 
12X,'DAMPING= 1 ,E12.6,/) 

DO 539 1=1 ,ITOT 
J=IPHI(I) 

IF(ITYPE(I) .EQ.2) GOTO 540 
PRINT 23,I,ZZ(J,K) 

GOTO 539 

540 PRINT 23,I,PHI(J,K} 

539 CONTINUE 

IF(NMEID.EQ.O) GOTO 548 
DO 47 1=1 , NMEID 
II=I+NME 

47 PRINT 23, IRESRE(I),ZZ(II,K) 

548 CONTINUE 
WRITE(2,19) 

PRINT 19 
PRINT 27 
WRITE(2,27) 

19 FORMAT (IX, 'DO YOU WANT TO REPEAT THIS WITH THE SAME FUNCTIONS? 
1' Y=1 N=0 ' ) 

READ(2,400)IT 
PRINT 301, IT 
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IF(IT.EQ.l) GOTO 888 
11 FORMAT( 2X , E12 , 6 , 2X , E12 . 6) 
END 


SUBROUTINE TRAN(F,IX,Y,NFR,J) 

COMMON/FREQ/ AM( 20) ,W(20) ,D(20) ,SM(20,20),IF0R(10),NF 
Y=0. 

DO 10 1=1, NFR 

Y=Y+( SM( J , I )*SM( IX , I )*2 .*D( I )*( F/W( I ) )/ 
1((1-(F/W(I))^2)**2+(2,*D(I) 1 * ( F/W(I))^2)) 

10 CONTINUE 
RETURN 
END 


SUBROUTINE FRF(ITW0,TW0M1,STFR,FIFR,NFR,I,AN0ISE,DSEED,KAY,IAVE) 
COMMON/ FREQ/AM( 20) ,W(20) ,D(20) ,SM(20,20) , IFOR( 10) ,NF 
C0MM0N/DATA/AI(256),AR(256),WW(256),F(10,512),UU(256) 

DOUBLE PRECISION DSEED 

* 

* DEVELOPES THE DISPLACEMENT FRFS 

* 

RANGE=FIFR-STFR 
DELT=RANGE/TW0M1 
NT I M= I TWO 
X=0.0 

DO 101 K=1 , ITWO 
KF0UR=K+(KAY-1)*ITW0 
KKK=KFOUR+I AVE* ITWO 
TK=FLOAT(K) 

WW(K)=STFR+(TK-1.)*DELT 

AI(K)=0. 

AR(K)=0. 

DO 101 KK=1,NF 
ISH=IFOR(KK) 

DO 101 J=1 ,NFR 
AK=AM(J)*W(J)*W(J) 

HI=(SM(I$H,J)*SM(I ,J)/AK)/ 

1((1.~(WW(K)/W(J))**2)**2+(2.*D(J)*WW(K)/W(J))**2) 

ANUM=(1.-(WW(K)/W(J))**2)*HI 

BNUM=-2.*D(J)*(WW(K)/W(J))*HI 

AI(K)=AI(K)+ANUM*F(ISH,KKK)+BNUM*F(ISH,KFOUR) 

AR( K)=AR( K)+ANUM* F( I SH , KFOUR)-BNUM*F( I SH , KKK) 

X=X+( AI( K)*AI( K)+AR( K)*AR( K) ) 

101 CONTINUE 


X=SQRT(X) 

C=AN0ISE*X/ITW0 
00 102 K=l, ITWO 
D$E=GGUBFS(DSEED) 

AN0I=2.*(2.*DSE-1. )*C 
AMAG=AI(K)*AI(K)+AR(K)*AR(K) 

AHAG=SqRT(AMAG) 

AI(K)=AI(K)+ANOI*AI(K)/AMAG 
AR(K)=AR( K)+ANOI *AR(K)/AMAG 
102 CONTINUE 
RETURN 

00 110 J=1,NTIM 
T=FLOAT(J) 

T=8.*ATAN(1.)*T/RANGE 
UU(0)=0,0 
00 110 K=1 , ITWO 
AN=FL0AT(K)-1, 

IF(K.EQ.l) AR(K)=AR(K)/2. 

UU(J)=UU(J)+AR(K)*C0S(T*DELT*AN)*2.-2.*AI(K)*SIN(T*DELT*AN) 
IF(K.EQ.l) AR(K)=2.*AR(K) 

110 CONTINUE 

WRITE(l.lll) (UU(J),J-1,NTIM) 

111 FORMAT! 4E14. 6) 

RETURN 

END 


SUBROUTINE MINV(AMAT,N,D,L,M) 
DIMENSION AMAT(900),L(29),M(29) 

* 

* MATRIX INVERSION SUBROUTINE 

"t 

42 CONTINUE 
D=1 . 

NK=-N 

DO 80 K=1 ,N 

NK=NK+N 

L(K)-K 

M(K)=K 

KK=NK+K 

BIGA=AMAT(KK) 

DO 20 J=K,N 
IZ=N*(J-1) 

DO 20 I=K,N 
IJ=IZ+I 

10 I F(ABS(BIGA)-ABS(AMAT( I J) ) ) 15,20,20 

15 BIGA=AMAT(IJ) 

L(K)=I 

M(K)=J 

20 CONTINUE 
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J=L(K) 

IF(J-K) 35,35,25 
25 KI=K-N 

DO 30 1=1, N 
KI=KI+N 

HOLD=-AMAT(KI) 

JI-KI+J-K 

AMAT(KI)=AMAT(JI) 

30 AMAT(JI)=HOLD 
35 I=M(K) 

IF(I-K) 45,45,38 
38 JP=N*(I-1) 

DO 40 J=1 ,N 

JK=NK+J 

JI=JP+J 

HOLD=-AMAT(JK) 

AMAT(JK)=AMAT(JI) 

40 AMAT(JI)=HOLD 

45 IF(BIGA) 48,46,48 

46 D=0.0 
RETURN 

48 DO 55 1=1, N 

IF(I-K) 50,55,50 
50 IK=NK+I 

AMAT(IK)=AMAT(IK)/(-BIGA) 

55 CONTINUE 
DO 65 1=1, N 
IK=NK+I 
HOLD=AMAT(IK) 

IJ=I-N 
DO 65 J=1,N 
IJ=IJ+N 

IF(I-K) 60,65,60 
60 IF(J-K) 62,65,62 
62 KJ=IJ-I+K 

AMAT( I J )=HOLD*AMAT( KJ )+AMAT( I J ) 
65 CONTINUE 
KJ=K-N 
DO 75 J=1,N 
KJ=KJ+N 

IF(J-K) 70,75,70 
70 AMAT( KJ )=AMAT ( KJ )/BIGA 

75 CONTINUE 
D=D*BIGA 

AMAT(KK)=1 ,/BIGA 
80 CONTINUE 
K=N 

100 K=(K-1) 

IF(K) 150,150,105 
105 I=L(K) 
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IF(I-K) 120,120,108 
108 0Q=N*(K-1) 
JR=N*(I-1) 

DO 110 J=1,N 
JK=JQ+J 
HOLD=AMAT(JK) 
JI=JR+J 

AMAT(JK)=-AMAT(JI) 
110 AMAT ( J I )=H0 LD 
120 J=M(K) 

IF(J-K) 100,100,125 
125 KI=K~N 

DO 130 1=1, N 
KI=KI+N 
HOLD=AMAT(KI) 
JI=KI-K+J 

AMAT(KI)=-AMAT(JI) 
130 AMAT(OI)=HOLD 
GOTO 100 
150 RETURN 
END 


SUBROUTINE AS I NE< FIST , I DTWO , STFR , FI FR , PH I , TW0M1 , XMEN ) 
DIMENSION XMEN(512) 

* 

* DEVELOPES TIME HISTORY OF SWEPT-SINE EXCITATION 

* 

DEL=FIST/TW0M1 
T0=1,/((TW0M1+1. )*DEL) 

DO 10 1=1 , I DTWO 
TI=FLOAT(I) 

F=STFR+DEL*(TI“1.) 

T=TO*( TI-1 . )/TWCMl 
10 XMEN( I )=( 3 . -T/TO)*SIN( F*T+PHI) 

RETURN 

END 


SUBROUTINE FFT(IDTWO,Z,TWOPI) 

COMPLEX Z(512) ,U,T,W 

* 

* DEVELOPES FRF OF THE SWEPT-SINE EXCITATION 

* 

11= I DTWO 
M=0 

PI=TW0PI/2. 

1 M=M+1 
11 = 11/2 


o tn 
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IF(II.GE.2) GOTO 1 
MV2=IDTW0/2 
NM1=IDTW0-1 
J=1 


DO 7 M,NM1 
IF(I.GE.J) GOTO 5 
T-Z(J) 

Z(J)“Z(I) 

Z(IK 

K*MV2 


IF(K,GEJ) GOTO 7 
J-J-K 


K=K/2 
GOTO 6 


J=J+K 

DO 20 L=1,M 
LE=2**L 
LEl=LE/2 
IH1.,0.) 

W=CMPLX(COS( PI/LEl ) »~SIN(PI/LE1)) 
DO 20 J==l » LEI 
DO 10 I=J,IDTWO,LE 
IP=I+LEl 


T=Z(IP)*U 

Z(IP)=Z(I)-T 

Z(I)»Z(I)+T 


u=u*w 


RETURN 

END 
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